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We carry out state-of-the-art optimization of a nuclear energy density of Skyrme type in the 
framework of the Hartree-Fock-Bogohubov (HFB) theory. The particle-hole and particle-particle 
channels are optimized simultaneously, and the experimental data set includes both spherical and 
deformed nuclei. The new model-based, derivative-free optimization algorithm used in this work 
has been found to be significantly better than standard optimization methods in terms of reliability, 
speed, accuracy, and precision. The resulting parameter set UNEDFpre results in good agreement 
with experimental masses, radii, and deformations and seems to be free of finite-size instabilities. 
An estimate of the reliability of the obtained parameterization is given, based on standard statistical 
methods. We discuss new physics insights offered by the advanced covariance analysis. 
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I. INTRODUCTION 

The goal of low-energy nuclear physics is to understand 
nuclei and how they react. This fascinating science prob- 
lem is relevant to other fields and to a gamut of societal 
applications. New vistas have been opened by experi- 
mental advances in the production of rare isotopes [l[ 
and new theoretical approaches [2| backed by unprece- 
dented computing power [3| . The rapid experimental de- 
velopments have resulted in a wealth of unique data from 
previously unexplored regions of the nuclear landscape. 
This situation poses a serious challenge to models of nu- 
clear structure and calls for their improved reliability and 
better-controlled extrapolability. 

Theorists seek to formulate a coherent framework for 
nuclear structure and reactions based on a well-founded 
microscopic theory that would deliver maximum predic- 
tive power with well-quantified uncertainties. To this 
end, the steady increase in computing power, currently 
crossing the petaflop barrier, has been beneficial. A 
paradigm for the new mode of nuclear theory is the Sci- 
DAC Universal Nuclear Energy Density Functional (UN- 
EDF) project [j|, an example of the close alignment of 
the physics research with the necessary applied mathe- 
matics and computer science research. 

This study is the fruit of such a partnership, under 
UNEDF. in which physicists collaborate with mathemati- 
cians and computer scientists on a specific science chal- 
lenge. Our long-term goal in UNEDF is to develop a 
spectroscopic-quality theoretical framework rooted in the 
nuclear density functional theory (DFT) Q. In the first 
phase of the project, we have developed efficient DFT 
solvers for the self-consistent Hartree-Fock-Bogoliubov 
(HFB) problem. Various improvements that we have im- 
plemented to carry out large-scale DFT calculations have 
been recently presented in [61, lZ[. These improvements 



enable comprehensive mass-table calculations, including 
all even-even nuclei and many different configurations in 
odd-even and odd-odd nuclei, in less than a day [^ Q . 

The second phase of the project concerns the devel- 
opment and optimization of the nuclear energy density 
functional (EDF). Since standard functionals are clearly 
too restrictive when one is aiming at a quantitative de- 
scription [13, [ill, the form of EDF needs to be im- 
proved. Novel functionals can be constructed from two- 
and three-nucleon interactions by using effective field the- 
ory and the density matrix expansion technique jl2Hl5l | 
and by using constraints from ab initio calculations for 
very light nuclei and nuclear matter. They can also be 
obtained by enriching density dependence and adding 
higher gradient terms in a systematic way J15l4l7l |. 

Having determined the form of the EDF, one must 
still optimize the coupling constants of the underlying 
energy density (ED). Indeed, all energy functionals, irre- 
spective of their theoretical foundations, rely on param- 
eters that must be directly fitted to experimental data. 
It has been realized recently that high-performance com- 
puting can positively impact the optimization strategy. 
Historically, most nuclear ED parameterizations, such as 
Skyrme or Gogny, were obtained by a direct fit to selected 
experimental data from finite nuclei and various nuclear 
matter properties (NMPs). Observables commonly in- 
cluded in the fit arc binding energies, proton radii, surface 
thickness, and/or single-particle (s.p.) energies of doubly 
closed-shell nuclei as well as NMPs (pseudo-observables) 
such as energy per particle of infinite and semi-infinite 
nuclear matter, saturation density, or incompressibility. 
This is the case, for example, for the SLy4 parameter- 
ization of the Skyrme functional of [18j . which we take 
(somewhat arbitrarily) as a reference point in our study. 
The Dl and DIS parameterizations of the Gogny interac- 
tion have also been obtained in such a framework (l9l.l20l. 



Wc refer to [5|, [2l|, ^^ for a more thorough discussion of 
various fitting strategies and protocols. 

In fact, very few examples of EDs are fitted to other 
types of data. For Skyrme EDFs only, we mention the 
early attempt of the SkM* parameterization [23] , which 
was adjusted semi-classically to account for the fission 
barrier of ^"^^Pu. The Brussels-Montreal set of EDFs 
has been optimized to data on deformed nuclei, although 
the actual fit is always performed with a spherical code 
by using a multistcp procedure. For example, in the 
early versions MSkl-MSk6, the deformation energy of 
the ground-state configuration was used to renormalize 
nuclear masses so that the optimization could proceed in 
spherical symmetry |2J] . Similarly, while in the later ver- 
sion HFB14, data on fission barriers were used, the core 
part of the fitting procedure was carried out in spherical 
geometry [2^. For SLy4 itself, several parameters were 
fixed at values empirically expected to yield a correct de- 
scription of giant resonance energy centroids in random- 
phase approximation (RPA) calculations, although no 
such calculation was included in the fit nor any quan- 
titative check performed a posteriori. 

The choice to restrict the data set of observables to 
those pertaining to nuclear matter and spherical nuclei 
has almost always been dictated by practical considera- 
tions: the cost of performing huge numbers of deformed 
HFB calculations was deemed too high. It was also 
rightly argued that the driving terms of the EDF could be 
pinned down by considering spherical nuclei only. With 
the need for more precision, however, the limitation to 
spherical nuclei and NMP is clearly not sufficient. The 
advent of supercomputers makes it possible to free our- 
selves from this restriction. 

Specifically, the availability of supercomputers has two 
consequences. First, one can now include in the set of fit 
observables data corresponding to deformed nuclei, odd- 
mass systems, excited states, and so forth. More com- 
prehensive data sets should better constrain the various 
channels of the energy functionals, for example, its de- 
formation or spin-polarization properties. It might soon 
become possible to directly optimize symmetry-restored 
EDFs [2g, either in a single-reference [27h30| or a mul- 
tireference j3lM33j framework. 

In addition, in our quest for improved EDFs, a key step 
is to understand various constraints imposed by experi- 
mental data on ED parameters and the resulting uncer- 
tainty margins. Early attempts to use statistical methods 
of linear-regression and error analysis j34| have been re- 
vived recently and applied to determine the correlations 
between ED parameters, parameter uncertainties, and 
the errors of calculated observables [10, [U [2l|, [22, [331 ■ 
This approach is essential for providing predictive capa- 
bility and cxtrapolability and for estimating the theoret- 
ical uncertainties. 

The purpose of this work is to revisit the problem of 
Skyrme ED optimization by (i) removing some of the pre- 
vious limitations with the help of modern computational 
resources, and (ii) applying regression diagnostics meth- 



ods on the resulting parameterization. To these ends, 
we perform functional optimization with a model-based 
method that is particularly adapted to costly function 
evaluations, such as when the objective function contains 
the result of hundreds of symmetry- unrestricted HFB cal- 
culations. In our model study, we focus on nuclear masses 
and radii, with a bias toward heavy nuclei. The final ED 
parameterization is subjected to a fully fiedged correla- 
tion and sensitivity analysis. While we do not claim to 
have found an end-all parameterization of the Skyrme 
EDF, we believe that the set of techniques we have ap- 
plied in this study can pave the way to a universal nuclear 
EDF of spectroscopic quality. 

The paper is organized as follows. In Sec.|ll]we briefly 
present the DFT framework used, in particular various 
parameterizations of the Skyrme EDF and their relations 
to nuclear matter properties. We also discuss the choice 
of experimental observables. Section lllll presents the spe- 
cific model-based algorithm used in this work and con- 
tains all the technical information related to large-scale 
HFB calculations. Results are discussed in Sec. [TV] Sec- 
tion |V| contains the conclusions of this work. 



II. THEORETICAL FRAMEWORK 

This section recalls the features of the Skyrme-DFT 
theory that are relevant to the optimization problem. A 
detailed presentation of the theory itself can be found 
in, for example, [5|, ISg . [37| and references therein. The 
main focus of the following discussion is on various pa- 
rameterizations of the Skyrme EDF and the selection of 
experimental observables chosen to constrain ED param- 
eters. 



A. Time-Even Skyrme Energy Density Functional 
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In nuclear DFT, the total energy of the nucleus is given 
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I 'H{r)(f'r 



(1) 



where % is the local energy density that is supposed to 
be a real, scalar, time-even, and isoscalar function of lo- 
cal densities and their derivatives. The Skyrme ED can 
be decomposed into the kinetic term, interaction ED x, 
pairing ED, Coulomb term, and additional corrections, 
such as the center-of-mass term. For the kinetic energy 
term, we set ;i2/2m=20.73553MeVfm2. The Coulomb 
Hartree term is calculated exactly, while the exchange 
term is computed by the Slater approximation. The con- 
tribution from the center-of-mass correction has the same 
structure as the kinetic term and leads to a renormaliza- 
tion of the nucleon mass \/m -^ (l/m)[l — 1/^]. All 
these prescriptions follow the SLy4 parameterization. 



The interaction ED ean be further decomposed into 
X = Xo + Xi, with 



+C^^''ptAp, 



Cf'ptV-Ju (2) 



where the isospin index t labels isoscalar (i=0) and 
isovector (t=l) densities. Since in this work we limit the 
discussion to even-even nuclei, the terms involving spin, 
spin-kinetic, and current densities [5|, 1371 138| are absent. 
The coupling constants Cf contain an additional depen- 
dence on the isoscalar density of the form 

CZ pI (3) 



cr 



CPP 



The standard Skyrme interaction ED therefore con- 
tains 13 independent parameters: 

{C«', CZ, Cl^', Cr. Cf , Cf '},=,,, and 7. (4) 

When dealing with the Skyrme interaction EDF (i.e., the 
functional that originates from the Skyrme interaction), 
the coupling constants ^ are uniquely related to the 
well-known (i, a;)-parameterization of the Skyrme inter- 
action 



{to, ii, i2, ^3, a;o, a;i, X2, ^3, to, io, &4, &4, 7}- 



(5) 



The equations connecting the C- and (i, x)- 
parameterization can be found, for example, in 

In this study, nucleonic superconductivity is described 
by the pairing ED: 



Xir) = E f 
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2 Po 



p\r), 



(6) 



where p is the local pairin g d ensity and pQ~0.16im ^ 
(mixed-pairing prescription |39|). 



B. Nuclear Matter Properties and Skyrme Energy 
Density Parameterizations 

The (t, x) and C-representations are natural param- 
eterizations of the Skyrme EDF, the former in terms of 
an effective, density-dependent two-body interaction and 
the latter as a general functional of the density. How- 
ever, these representations do not provide a straightfor- 
ward connection to physical observables; hence, it is not 
immediately obvious what the search range for these pa- 
rameters should be. It is therefore advantageous to relate 
them to fundamental properties of symmetric and asym- 
metric homogeneous nuclear matter, which have a clear 
tsical interpretation and the range of which is known 
MM- 

The starting point in the discussion of NMPs is the 
equation of state (EOS) of the infinite homogeneous nu- 
clear matter: E/A = W{pn, Pp)- The Coulomb energy is 



disregarded, all gradient terms vanish, and the kinetic en- 
ergy density is replaced by its Thomas-Fermi expression. 
Assuming an unpolarized system, one can also ignore 
terms involving time-odd spin densities and currents. 

The expansion of W{pn,Pp) around the equilibrium 
density pc and / = can be written as 

Wipn, Pp) = W{po, I) = Wipo) + S2{po)I^ + 0{I^), (7) 

where / = pi/ po = [pn — Pp)/ Po is the relative neutron 
excess, po = Pn + Pp, pi = Pn - Pp, 



'NM 



W{po) = — ;— + ^^ (Po - Pc) + -— ^ (po - Pc)^ , (8) 
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(9) 



In these equations, E^^/A stands for the total energy 
per nucleon at equilibrium, p^M represents the nucle- 
onic pressure, K^^ is the nuclear matter incompressibil- 
ity, af^ is the symmetry energy coefficient, Lf^ rep- 
resents the density dependence of the symmetry energy, 
and AK^^ is a correction to the incomprcssibility. 



1. Symmetric nuclear matter 

In the regime of symmetric nuclear matter (SNM), 
Pn = Pp = Po/2 and /=0, which eliminates all isovec- 
tor terms. The isoscalar kinetic energy density is 

.2 \ 2/3 



To = CkPo , 



a = 



3 /37r^ 



(10) 



The nuclear matter saturation curve M^(po) is expected 
to have the following properties: 



pNM 
A 
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0.16 fm^^ 
2dW{po) 



dpo 
W{pc) 



0, 



(11) 
(12) 



PO=Pc 



- 16 MeV. (13) 



The value of the incompressibility modulus is related to 
the centroid energies of giant isoscalar monopole reso- 
nances in isospin-symmetric nuclei ^^ and is expected 
to be El H 



.NM_ a^2d^W{po) 



if"" = 9p^- 



dPo 



220 ± 10 MeV, 



(14) 



Pa=Pc 



with a strong preference for 230 MeV (45|. Another im- 
portant NMP, entering the SNM EOS indirectly, is the 
isoscalar effective mass 



AC = 



2m dE 
h? dTQ 



(15) 



PO=Pc 



which quantifies the momentum-dependence of the mean 
field and drives the density of the s.p. spectrum. An 
appropriate value for a fit to experimental s.p. energies is 
M* = 1 [4a] , while ab initio calculations performed at the 
Brueckner-Hartree-Fock level in INM suggest a slightly 
lower value for the Landau (Fermi-level) effective mass 
extracted from the on-shell s.p. spectrum |47H5(]| |. Mass 
fits also seem to favor a value close to unity, although 
significant freedom exists [51[ . 

The SNM EOS expressed in terms of the coupling con- 
stants of the Skyrme EDF is 
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2/3 



(16) 



Computing the quantities (|13m5p using ((TB)) allows us 
to express the coupling constants Cqq , Cq^, Cq^ and the 
power 7 in terms of E^^/A, P^m ^ q^ ^nm ^^^ M*-\ 
The resulting expressions are [40] as follows: 



GSS = ^{l^[(2-37)Mr^-3] 
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K 



NM 



£,{6Mr^^9)r. + 9^ 



(17) 
(18) 

(19) 
(20) 



where r^ 



Ckpi 



5/3 



The value of S2{pc) varies from 28 to 36 MeV among 
EDFs extrapolated to nuclear matter [4l|, [52| . It is un- 
derstood [53l l54| that nuclear masses constrain a com- 
bination of the symmetry- and surface-symmetry energy 
parameters in a given EDF, and this fact explains the 
large spread of values. 

The variation of the density-dependent symmetry en- 
ergy 52 with po is usually parameterized through 

dS2{po) 



'NM 



T- IN iVi o 



dpo 



(26) 



the value of which appears correlated with the thickness 
of neutron skins in asymmetric nuclei (see |35| and ref- 
erences therein) . An empirical determination of this pa- 
rameter yields L^yM ^ 80 ± 30 MeV [H, [Hj. One now 
introduces 

d^S2ipo) 



AK'''^' = 9pi 



d'^po 



(27) 



which affects the incompressibility of the ANM and thus 
the isoscalar monopole resonance energies in neutron-rich 
nuclei ^. For the SLy4 EDF, the values of the last 
two parameters were determined by the fit to the neu- 
tron matter EOS. We let these quantities be constrained 
by our experimental data set. We will see whether 
these data leave enough freedom to apply additional con- 
straints in the regression analysis. 

The momentum dependence of the mean field is also 
affected by isospin: neutron and proton effective masses 
are different in asymmetric matter JSTJ , an effect quanti- 
fied by the isovector effective mass 

, ^* 1 , r* 1 ^m dE 

M = M 

ft2 d^^ 



Po 
I = 



Pc 



(28) 







The EOS of homogeneous asymmetric nuclear matter 
can be written as 



2. Asymmetric nuclear matter 

In asymmetric nuclear matter (ANM), neutron and 
proton densities are different, and isovector terms are 
nonzero. The local and kinetic energy densities are 



Pi 
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F±il) 
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CypTF+{I), 
C^pl'^F^il), 
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(21) 
(22) 
(23) 
(24) 



The nuclear matter EOS W{I,p) now depends on the 
relative neutron excess /. The most important parameter 
characterizing the isospin dependence of the ANM EOS 
is the symmetry energy at saturation density. 



^. , NM_ l d^W{p,J) 

^2\Pc) — flsym — 2 'Jp 
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(25) 



+crckP^/'/F_(/) 

+ [CZ + CZpI + I' (Cfo' + C^dP'o)] po. (29) 
Just as for SNM, we compute the quantities (|25ll28p from 
pg]) and obtain an expression for Cfp , Cf^ , and Cf^ [401 ■ 
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(30) 



(31) 



(32) 



Using relations (|20|) - ((32)) . we express 7 of the original 
13 parameters Qi of the Skyrme EDF as functions of nu- 
clear matter properties. The remaining 6 are not known 
exactly and should therefore not be used as rigid con- 
straints [22|. However, the expected values of all these 
NMPs are sufficient to provide well-defined intervals of 
variation during the optimization process. The 6 remain- 
ing coupling constants arc the isoscalar and isovector 
Cf '' , spin-orbit Cf ' , and tensor C/ terms. Conse- 
quently, the Skyrme EDF depends on the following 13 
parameters: 



{ 



„ pNM/4 71 r* r<-NM „NM 7-NM 



M:,C^^Cf^^C^^Cf^^Co^^cf } . (33) 



C. Fit Observables 

To calibrate the EDF, we selected a pool of fit observ- 
ables that constitute the UNEDF experimental database 
J58| . The purpose of the database is to provide a stan- 
dard and comprehensive set of experimental data that 
can be used to systematically optimize EDFs. Since 
we wish to provide, together with the optimized set of 
parameters, a measure of its intrinsic quality via the 
error and sensitivity analysis, for every observable an 
error bar should also be defined. We organized our 
database into three major categories - spherical, de- 
formed, and symmetry-unrestricted - which reflect the 
level of symmetry-breaking of the underlying EDF and 
thereby the complexity of its numerical implementation. 
More details can be found in |58j . 

The focus of this work is on a well-controlled optimiza- 
tion methodology, and the emphasis is on global nuclear 
properties such as masses and proton radii. Our func- 
tional is therefore restricted to time-even densities, and 
only spherical or axially deformed nuclei are considered. 
The chosen observables embrace data for 72 nuclei, which 
are proven to allow a reasonable DFT description. The 
selected experimental data set is presented in Fig. [T] As 
can be seen, the emphasis is on the heavy nuclei. In- 
deed, there are only 11 nuclei with A<66 in our data 
set. Below, we give a detailed description of the set of fit 
observables used in this work. 
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FIG. 1: (color online) Experimental set of fit observables used 
in this work. The set contains data for 11 nuclei with A<66 
and 61 nuclei with A>106. 



TABLE I: Nuclear binding energies (in MeV; the electronic 
energy correction has been subtracted) [53| for the 44 de- 
formed nuclei selected in this work. The column marked "#" 
is the data point number. 



# 


Z 


N 


E 


# 


Z 


N 


E 


1 


108 


156 


-1925.697 


23 


94 


144 


-1800.523 


2 


106 


154 


-1908.038 


24 


92 


144 


-1789.701 


3 


104 


152 


-1889.709 


25 


92 


142 


-1777.858 


4 


102 


154 


-1897.729 


26 


90 


142 


-1766.015 


5 


102 


152 


-1884.685 


27 


72 


104 


-1418.407 


6 


102 


150 


-1870.386 


28 


70 


108 


-1431.260 


7 


100 


156 


-1901.673 


29 


70 


100 


-1377.760 


8 


100 


154 


-1890.112 


30 


68 


104 


-1391.213 


9 


100 


152 


-1878.056 


31 


68 


102 


-1378.695 


10 


100 


150 


-1864.657 


32 


66 


102 


-1362.591 


11 


100 


148 


-1850.682 


33 


66 


100 


-1350.474 


12 


100 


146 


-1836.305 


34 


66 


98 


-1337.714 


13 


98 


156 


-1891.281 


35 


66 


96 


-1323.785 


14 


98 


154 


-1880.445 


36 


66 


94 


-1309.134 


15 


98 


152 


-1869.165 


37 


66 


92 


-1293.725 


16 


98 


150 


-1856.954 


38 


66 


90 


-1277.701 


17 


98 


148 


-1843.959 


39 


64 


98 


-1321.473 


18 


98 


146 


-1830.429 


40 


64 


96 


-1308.992 


19 


98 


144 


-1816.428 


41 


64 


94 


-1295.597 


20 


96 


150 


-1847.037 


42 


64 


92 


-1281.300 


21 


96 


148 


-1835.059 


43 


64 


90 


-1266.329 


22 


96 


144 


-1809.502 


44 


64 


88 


-1251.187 



1. Deformed nuclei 



In our optimization, we considered binding energies 
of 44 well-deformed even-even nuclei shown in Tabic ID 
Candidates were selected from an HFB mass-table cal- 
culation with the SLy4 parameterization requiring that 
their ground-state equilibrium deformation be greater 
than |/3[=0.25. Since the majority of atomic nuclei are 
deformed in their ground states, by including binding en- 
ergies of deformed systems in the database, one hopes to 
better probe the surface properties of the EDF. 



2. Spherical nuclei 

Table HIl lists the nuclear masses of a selected set of 28 
spherical nuclei considered in the fit. In these nuclei, cor- 
relations beyond mean-field are expected to be relatively 
constant [23. Since the list includes doubly magic nuclei. 



it should provide strong constraints, as these nuclei tend 
to deviate from global mass trends [1^. Moreover, the 
masses of ^°Ca, ^^Ca, and ^^Ni help constrain the spin- 



orbit term [57|, |60|, |6l| . All the masses of spherical and 
deformed nuclei given in Tables Ulllll have been corrected 
for the electronic binding energy. The nuclear binding 
energy Enuc{Z, N) is given by 



E,-,uc{Z,N)=E,to{Z,N)-E, 



(34) 



where Es,to{Z, N) is the atomic binding energy and E^i 
-1.433 X 10-5^2-39 MeV. 



TABLE II: Nuclear binding energies (in MeV; the electronic 
energy correction has been subtracted) [62| for 28 spherical 
nuclei selected in this work. The column marked "#" is the 
data point number. 



# 



N 



E 



# 



N 



E 



45 
46 
47 
48 
49 
50 
51 
52 
53 
54 
55 
56 
57 
58 



82 
82 
82 
82 
82 
82 
82 
82 
82 
50 
50 
50 
50 
50 



132 
130 
128 
126 
124 
122 
120 
118 
116 
74 
72 
70 



-1662.762 
-1653.988 
-1645.030 
-1635.909 
-1621.803 
-1606.984 
-1591.666 
-1575.833 
-1559.483 
-1049.835 
-1035.365 
-1020.375 
-1004.785 
-988.535 



59 
60 
61 
62 
63 
64 
65 
66 
67 
68 
69 
70 
71 
72 



50 
50 
50 
28 
28 
28 
28 
28 
20 
20 
20 
20 
20 
20 



64 
62 
58 
36 
34 
32 
30 
28 
30 
28 
26 
24 
22 
20 



-971.406 
-953.335 
-914.424 
-561.714 
-545.217 
-526.801 
-506.459 
-483.949 
-427.473 
-415.972 
-398.751 
-380.942 
-361.877 
-342.033 



For the same 28 spherical nuclei, we also consider the 
proton rms point radius (i?p), which we extract from the 
charge radius (i?ch) of l£3 using the standard relation: 



(i?,\) = (i?2) + (r2) + 9(r2), 



N 

'z' 



(35) 



where the proton charge radius, 
the neutron charge radius, (r^ 




= 0.877 fm, and 
0.1161 fm^, were 



taken from [64j . The values of proton radii used in this 
work are listed in Table IIIII 



Since the particle-hole and particle-particle channels 
cannot easily be disentangled, we must also include ob- 
servables that will help us pin down the magnitude of 
pairing correlations. Usually, the pairing part of the 
EDF is constrained by considering the odd-even stagger- 
ing (OES) of binding energy (see [6a] for a recent survey). 
Additional constraints on the pairing ED may be imposed 
by taking calculated pairin g ga ps in symmetric nuclear 
matter and neutron matter [66| . This strategy has been 
adopted by the Brussels-Montreal group in their most 
recent model HFB-17 [eTJ . 

In this work, we constrain pairing EDF by means of 
the OES defined by a 3-point formula A^^) [el, [el]. As 



TABLE III: Proton rms radii (in fm) [631 for the 28 spherical 
nuclei selected in this work. The column marked "#" is the 
data point number. 



±_ 
73 
74 
75 
76 
77 
78 
79 
80 
81 
82 
83 
84 
85 



N 



# 



N 



82 
82 
82 
82 
82 
82 
82 
82 
82 
50 
50 
50 
50 
50 



132 
130 
128 
126 
124 
122 
120 
118 
116 
74 
72 
70 
68 
66 



5.506 
5.488 
5.469 
5.450 
5.439 
5.428 
5.418 
5.403 
5.394 
4.609 
4.598 
4.586 
4.573 
4.558 



87 
88 
89 
90 
91 
92 
93 
94 
95 
96 
97 
98 
99 
100 



50 
50 
50 
28 
28 
28 
28 
28 
20 
20 
20 
20 
20 
20 



64 
62 
58 
36 
34 
32 
30 
28 
30 
28 
26 
24 
22 
20 



4.542 
4.527 
4.492 
3.787 
3.765 
3.733 
3.689 
3.661 
3.437 
3.390 
3.412 
3.432 
3.420 
3.382 



TABLE IV: Values of the neutron and proton average odd- 
even mass staggering (in MeV) considered in this work. The 
column marked "#" is the data point number. 





Neutrons 






Protons 




# 


Z 


N 


X(3) 


# 


Z 


N 


X(3) 


101 


100 


152 


0.515 


105 


96 


148 


0.566 


102 


92 


144 


0.569 


106 


92 


142 


0.606 


103 


72 


104 


0.675 


107 


68 


102 


0.504 


104 


66 


98 


0.679 


108 


66 


94 


0.728 



customary, the theoretical result for even particle num- 
ber TV is compared with the experimental A^"^) for iV -I- 1 
[63 | . We took four values of A'^'^^ for neutrons and four for 
protons; sec Table HVl All these nuclei belong to the de- 
formed set of Table HI Our choice has been motivated by 
the observation that fitting pairing properties in spheri- 
cal systems, where the level density is much greater, may 
lead to an underestimation of the overall pairing strength 

With the fairly simple pairing ED ^ that we use, it is 
not essential to require very high precision for the OES. 
For that reason, in order to be free from local fluctu- 
ations, we chose in each even-even nucleus the average 
over the two even-odd or odd-even isotopes: A„ ' (N) = 
[A(3)(7V-l) + A(^^)(A^+l)]/2. Including average values of 
A*^^' in our data set ensures that the magnitude of pair- 
ing correlations is correct and remains such throughout 
the fitting procedure. Theoretical OES values have been 
computed from the average HFB pairing gap [70, IZll • 



III. OPTIMIZATION ALGORITHM 

This section briefly presents the new algorithm used 
in our optimization. We refer to it by the acronym 



POUNDerS, standing for Practical Optimization Using No 
Derivatives (for Squares) . We also provide the numerical 
parameters used in the HFB calculations, and we give 
the characteristics of the objective function used in the 
optimization. 



A. Derivative-Free Optimization Method 

To outline our algorithm, we adopt the following nota- 
tion. We denote the set of parameters/coupling constants 
of the Skyrme EDF to be fitted by x G R"" , where rix is 
the number of coupling constants of components Xk to fit. 
We define a composite fit function made of Dt different 
types of data: nuclear masses, proton radii, and so on. 
The number Ui of data points for a given type i may vary; 
for example, we have more masses than rms radii. The 
output of the calculation for type i is denoted by Sij(x) 
for nucleus j and obviously depends on the parameteri- 
zation of the functional, that is, the vector x G M"^ . For 
type i and nucleus j, the experimental value of a given 
observable is denoted dij. 

While many objectives are possible, we minimize the 
weighted sum of squared errors 



X'(x) 






Ud -n. 






i=i j=i 



(36) 



where n^ = X^iJ^i "-* denotes the total number of data 
points being fit. The weights Wi > render the type i 
difference dimensionless and are chosen to balance the 
goals of fitting different observable types simultaneously. 
The objective ([55)) is a special case of the nonlinear 
least squares function 



1 "'" 1 

/(x) = -5]j^.(x)2==-||F(x)| 



(37) 



where the function F : M"" -^ M""* yields the vector of 
reduced errors. Most optimization approaches to mini- 
mizing ([57]) are based on Newton's method, whereby / is 
replaced by its second-order expansion 



numerical differences. In such a case, the optimization 
algorithm must be derivative-free, relying only on the 
function value outputs F(x). Popular algorithms in this 
setting include the Nelder-Mead (N-M) method and other 
direct search algorithms [T^l and genetic algorithms and 
other heuristics (73|. However, a recent benchmarking 
study |74[ found that methods that form a smooth ap- 
proximation model of the objective in order to exploit the 
smoothness and structure of the objective may be able 
to obtain better solutions in fewer evaluations. 

In the case of nonlinear least squares, we follow the 
approach of forming a quadratic model for each compo- 
nent. 



qii-^ + S) 



^,(x)+^'^g, + ^J^H,<5, 



(39) 



with gi and H^ = H^ playing the role of the unknown 
derivatives Vi^i(x) and V^Fj;(x), respectively. We ob- 
tain the model parameters gi and H^ by requiring that 
the model g^ agree with the true function Fi on a set 
<Y of X values at which Fi is known. Mathematically, 
these parameters are solutions to the convex quadratic 
program 



min {||H,||j. : g,(xfe) = F,;(xfc) Vx^ G X} , 



(40) 



where || \\f is the Frobenius norm and the interpolation 
set X contains between n^ + 1 and {ux + l){nx + 2)/ 2 
points satisfying geometric conditions detailed in (75|,|76|. 

The quadratic model qi cannot be expected to approx- 
imate Fi at X values far from the points in X. Hence, 
we use a trust region framework, whereby the model qi 
is trusted only close to a base-point x. Given a radius 
A > 0, we let B = {x e K"" : ||x - xj| < A} denote the 
spherical neighborhood within which we trust qi. Corre- 
spondingly, the interpolation points in X should not be 
too far away from B. 

Provided that we know the entire vector of observables 
F(xfe) at each Xfc G X, we can obtain a set of model 
parameters {(gi,Hi)}"^^, which we trust inside a com- 
mon region B centered about x. We can thus form a 
derivative-free model of the quadratic 



/(x + ^)«/(x) + J^J(x)^F(x) 

+ i^^ ( J(x)^ J(x) + £ F,;(x)V2f,(x) J S, (38) 

where S G M"°= and J(x) is the Jacobian matrix J(x) = 
[V^i(x),-.-,V^„,(x)F. 

In the problem at hand (and many others), the deriva- 
tives of Fi(x) with respect to x, Vi^i(x) and V^Fi(x), 
exist for virtually all x, but their calculation for use in 
the optimization is impractical. Indeed, although deriva- 
tives of binding energies can be obtained through the 
Feynman-Hellman theorem, other observables such as 
radii would require the use of perturbation theory or 
a cumbersome and potentially imprecise calculation of 



m{^ + S) = f{i) + S^Y.F,{±)g, 

i=\ 
1 "rf 

+ -<5^^(g,gf + F,(x)H,)J. (41) 

Since we trust this model within B, we expect that a 
better x can be obtained by solving the trust region sub- 
problem min5{m(x-|-^) : x+ J G B}. This problem mini- 
mizes a quadratic with known derivatives over a compact, 
convex region and is hence decidedly easier than the orig- 
inal problem. The observables are then evaluated at the 
solution to this subproblem so that we obtain F(x -|- 5). 
An iterative Newton-like procedure is thus obtained. 
We note that the trust region radius A grows and shrinks 



from one iteration to the next depending on the ratio of 
the actual decrease obtained at the new point versus the 
decrease predicted by the model in (|4T|) . Similarly, our 
current estimate of the solution, x, is changed only if an 
adequate decrease of the function was obtained or if we 
achieved a simple decrease in the function value and the 
geometry of the interpolation set X gives us confidence. 
If we did not adequately decrease /, we must evaluate at 
an additional x value in order to improve the geometry 
of the set X in subsequent iterations. 



B. Numerical Parameters 

The evaluation of the function ([M)) at point x requires 
72 HFB calculations to generate the Sij{x.) points for 
the 72 nuclei j taken in the data set. All HFB calcula- 
tions were performed with the code hfbtho [73] ■ This 
code solves the Skyrme-HFB equations in the harmonic 
oscillator (ho) basis assuming axial and reflection sym- 
metry. In our optimization, we used a spherical basis of 
^^811011=20. The oscillator frequency was determined for a 
given nucleus of mass number A according to the formula 
/iWoscii = 1.2 X -^jj MeV [23 • These two choices guaran- 
tee good convergence of the HFB energy with respect to 
the basis size, within about 150 keV of the exact value 

Pairing correlations were described by the pairing ED 
^ with different pairing strengths for protons and neu- 
trons, VfP 7^ Vq . As customary for zero-range pairing 
forces, a cut-off of i5™^60MeV is used to truncate the 
quasi-particle space [70]. In order to avoid pairing col- 
lapse, the Lipkin-Nogami prescription was systematically 
applied according to [80{ . 

Taking into account the 13 parameters of the Skyrme 
EDF and the 2 additional parameters in the pairing chan- 
nel requires a 15-parameter search. We have made two 
additional simplifications. First, the tensor coupling con- 
stants Cq and C/" were set to 0. This choice was moti- 
vated by our requirement to take as a reference point the 
original SLy4 parameterization of [l8| where these terms 
were not included. Second, preliminary tests indicated 
that the isovector effective mass was poorly constrained 
by our data set. As a result, the obtained values of M* 
were clearly nonphysical with regard to the discussion in 
[Slf . In the final run we therefore discarded M* from the 
list of free parameters and kept the original SLy4 value. 

The final optimization was therefore carried on a set 
of 12 parameters (10 for the Skyrme ED plus 2 pairing 
strengths): 



/„ pNM//i T^NM „NM rNM 71, r 



*-l 



c, 



pAp ^pAp^ yr.^ yp^ ^pVJ^ ^pVJ^ ^ (42) 



with Ci^ = Cf = and M*-i = 1.249. 

For scaling purposes, the optimization algorithms 
tested here require the domain of variation of the var- 
ious parameters x to be specified. Since, in practice, a 



large subset of x represents symmetric and asymmetric 
nuclear matter properties, the range of variation can be 
easily set up, even if the exact values are not known. Ta- 
ble |V| in Sec. IIV All lists the scaling intervals adopted in 
our optimization. 

Following the discussion in Sec. IlICi our objective 
function p6p contains Z?t= 3 data types: nuclear masses 
(i=l), proton rms radii (i=2), and OES differences («=3). 
The total number of data points is 71^=108 and breaks 
down into ni=72 nuclear masses (28 spherical and 44 
deformed), ri2=28 rms proton radii, and n3=8 OES dif- 
ferences (4 for neutrons and 4 for protons). The values 
dij of the experimental data points are given in Sec. Ill Cl 

The weights Wi in the objective function are used to 
render all quantities dimensionless and to allow for a com- 
posite y^ function. The weights must be chosen so that 
all reduced errors are of the same order of magnitude: 
they reflect the expected theoretical uncertainty that one 
can assign to a given observable, which is generally larger 
than the corresponding experimental uncertainty for our 
data set. In the optimization described here, we chose 
Wmass=2.0MeV, uiradii=0.02fm and ii;oES=50keV. 



IV. RESULTS 

This section contains the optimization results. In Sec. 
IIV Al various properties of the resulting ED parameter- 
izations are explored. Section flV Bl illustrates the versa- 
tility of our approach by providing a detailed correlation 
and sensitivity analysis. 



A. Optimized Functionals unedfnb and unedfpre: 
Properties and Stability 

We give in this section the final parameterization of the 
Skyrme functionals UNEDFnb and UNEDFprc that mini- 
mizes the x^ objective function p6p . and we perform a 
number of checks to probe the quality of the resulting 
functionals. In particular, we test their stability with the 
RPA response function, check that both the spherical 
and deformed shell structure are on par with other pa- 
rameterizations, and discuss various global performance 
indicators. 



1. Solution to the optimization problem 

The optimization of a nuclear energy functional is a 
complex problem. The objective function is the com- 
pound result of many different full HFB calculations, 
each the result of a self-consistent iterative procedure. 
In principle, such a function lends itself naturally to par- 
allelization, although the different times of calculation of 
spherical and deformed configurations requires fine load 
balancing. Overall, the cost of one function evaluation 
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FIG. 2: (color online) Convergence of the minimization of Eq. 
(|36|l with the standard Nelder-Mead algorithm (dashed) and 
the model-based POUNDerS (solid line). 



TABLE V: Values x of the optimization parameters x of Eq. 
(|42p at the solution with no bounds imposed (Skyrme func- 
tional UNEDFnb). pc is in fm^^; E^^'^/A, K^^ , af^, and 
Lfy^ are in MeV; 1/M* is dimensionless; C^^" and C^'' in 
MeVfm^; and Vq and VJf in MeVfm^. The range of variation 
provided to the optimization is shown in the column "Scaling 
Interval," the initial values in column x(""*)^ and the final 
values in x'*""''. 



k 


X 


Scaling Interval 


x(init.) 


X(fln-) 


1. 


Pc 


[ -f0.14 , +0.18] 


+0.160 


0.151046 


2. 


^NM/^ 


[ -17.00, -15.00] 


-15.972 


-16.0632 


3. 


iv-NM 


[+170.00, +270.00] 


+229.901 


337.878 


4. 


NM 


[ +27.00, +37.00] 


+32.004 


32.455 


5. 


rNM 
-t^sym 


[ +30.00, +70.00] 


+45.962 


70.2185 


6. 


1/m: 


[ +0.80, +2.00] 


+1.439 


0.95728 


7. 


QpAp 


[-100.00, -40.00] 


-76.996 


-49.5135 


8. 


QpAp 


[-100.00, +100.00] 


+15.657 


33.5289 


9. 


Vo" 


[-350.00, -150.00] 


-258.200 


-176.796 


10. 


Vo' 


[-350.00, -150.00] 


-258.200 


-203.255 


11. 


c^^-' 


[-120.00, -50.00] 


-92.250 


-78.4564 


12. 


cf' 


[-100.00, +50.00] 


-30.750 


63.9931 



can typically amount to 10 minutes on a standard com- 
puter cluster. With such costly evaluations, the number 
of evaluations required to minimize p6p can rapidly be- 
come an issue. 

In addition, we have no prior knowledge of the multidi- 
mensional surface of the objective function in the param- 
eter space. There is no guarantee that the parameters are 
all independent, and, as we show later, there are correla- 
tions between them, which make the topography of the 
surface complex. 

These observations suggest that two important fea- 
tures of a good optimization algorithm should be the 
speed of convergence and the ability to converge to a 
true minimum, if only a local one, without being mis- 
led by narrow valleys and saddle points. Figure [5] shows 
the performance of the standard Ncldcr-Mead (N-M) al- 
gorithm, as implemented in the TAO code |82| . on our 
objective function, compared with the new model-based 
algorithm presented in Sec. IIII Al We note that the 
POUNDerS method attains a value of y^ close to the final 
one after only 25-30 iterations, whereas after more than 
300 iterations the N-M algorithm yields a solution that 
is still a factor of 2 away. Moreover, there seems to be 
a stagnation of the N-M method at around 15-65 itera- 
tions, which may prematurely suggest that the minimum 
has been found. Yet, in this plateau the x^ is still about 
5 times larger than at the final solution. 

Table |V] shows the values of the optimization parame- 
ters ((1^ at the solution (dubbed UNEDFnb in the follow- 
ing). The starting values were given by the SLy4 param- 
eterization. The most notable change affects the effective 
mass: starting from M* « 0.7, the final value is close to 
1 . which ensures a level density more compatible with the 
empirical one (even though there is no obvious reason for 
this to happen, given the data set employed). As will be 



discussed in Sec. IIV A 61 without being steered, the opti- 
mization gives the correct hierarchy of pairing strengths, 
namely, \Vq\ > \Vq\, to reflect the missing momentum- 
dependence and Coulomb contribution, as pointed out in 

iiiiilll- 

A standard measure of the quality of the optimization 
is the rms deviation (RMSD) of data type i at the solu- 
tion x: 



RMSD(i) = 



\ 



1 "' 



(43) 



j=i 



For our set of fit observables, the RMSDs for 
various types of data are RMSD(mass)=0.966MeV, 
RMSD(radu)=0.014fm, and RMSD(OES)=57keV. For 
comparison, the value of RMSD(mass) for SLy4 on the 
same data set is 9.95 MeV. 

A close examination of Table |V] shows that, while 
most of the parameters of UNEDFnb have values in 
the normally accepted range, the incompressibility 
iir^^=338MeV is far too large. This would seriously 
limit the usability of UNEDFnb in nuclear structure cal- 
culations, in particular in studies of collective modes such 
as monopole vibrations. 

We therefore performed another minimization, using 
the same scaling intervals, but imposing hard bounds on 
the NMPs. A similar strategy was adopted in [8g, where 
hard bounds on K^^ were imposed during optimization 
of BSkl3 EDF. 

Tabic I VII displays the parameter values of the Skyrme 
functional UNEDFpre optimized in such a way. At con- 
vergence, the nuclear incompressibility and scalar effec- 
tive mass appear at their respective bounds of 230 MeV 
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TABLE VI: Same as Table |V] but for the case with bounds 
(Skyrme functional UNEDFpre). 



k 


X 


Bounds 


x(init.) 


X(fln.) 


1. 


Pc 


[+0.15,+0.17] 


+0.160 


0.160526 


2. 


^NM/^ 


[-16.2,-15.8] 


-15.972 


-16.0559 


3. 


^^NM 


[-hl90, -h230] 


+229.901 


230 


4. 


NM 
"'sym 


[ +28, +36] 


+32.004 


30.5429 


5. 


rNM 
-^syin 


[ +40, +100] 


+45.962 


45.0804 


6. 


1/m; 


[ +0.9, +1.5] 


+1.439 


0.9 


7. 


^pAp 


[— O0,+CX)] 


-76.996 


-55.2606 


8. 


QtpAp 


[—00, +CX)] 


+15.657 


-55.6226 


9. 


Vo" 


[— OQ, +CX)] 


-258.200 


-170.374 


10. 


vi 


[— CXD, +0o] 


-258.200 


-199.202 


11. 


c'S^' 


[— oo, +oo] 


-92.250 


-79.5308 


12. 


^PVJ 


[— oo, +oo] 


-30.750 


45.6302 



and 1.11 (1/M* = 0.9); that is, these NMPs are ac- 
tively constrained. The rms deviations obtained for UN- 
EDFpre on our set of fit observables are still respectable: 
RMSD(mass)=1.455MeV, RMSD(radii)=0.016fni, and 
RMSD(OES)=59keV. 



2. Stability check of UNEDF716 and UNEDFpre 

It is known that some Skyrme ED parameterizations 
are prone to finite-size instabilities [8l|, l86l - l88j . For in- 
stance, in the time-even channel, the term Ci ''piApi 
can lead to divergences of the HFB iterative procedure. 
When searching for new functionals, it is therefore cru- 
cial to test comprehensively the stability of the functional 
parameterization. Here, the RPA linear response theory 
[89| is the tool of choice. The full RPA response in infi- 
nite matter has been derived for Skyrme EDFs [90l - l92l l. 
and applications pertaining to the stability of Skyrme 
functionals have been reported in [3, |8l[ . 

Without entering into details, a general expression for 
the RPA response function n(w, q) in SNM can be writ- 
ten as ii 



n(w,q) 



4no(a;,q) 



(44) 



where w is the excitation energy; q is the transferred 
momentum (or wave number of the density fluctuation) ; 
no(a;,q) is the noninteracting response (or Lindhard 
function); and D[uj,q^ is the dielectric function, equal 
to unity in noninteracting SNM. 

The value n(w = 0, q) corresponds to the static suscep- 
tibility of the system to finite-size perturbations. With 
the above sign convention, n(w = 0, q) should be pos- 
itive for all values of q and the density po- A change 
of sign with either variable corresponds to £'(w, q) = 0; 
hence, the occurrence of a pole indicating the existence of 
a zero-energy collective mode. In the isospin channel, the 
short-wavelength (high-q) behavior is driven essentially 




FIG. 3: (color online) Dielectric function D{u = 0, q) for the 
scalar-isovector channel in SLy4, UNEDFpre and UNEDFnb as a 
function of the transferred momentum g, for kp ~ 1.33 fni~^. 



by the combination of coefficients Cf — Cf ''q^ |57l.[90|. 
The magnitude of the latter correlates well with the oc- 
currence of instabilities in calculations of finite nuclei. 

Figure 13] shows the dielectric function D{uj = 0, q) as a 
function of 5, in the scalar-isovector perturbation chan- 
nel, at saturation density in SNM. When D{uj = 0, q)=0, 
finite-size instabilities could potentially develop and hin- 
der the usability of the functional. This situation does 
not occur for UNEDFnb and UNEDFpre, which yield a di- 
electric function even more "stable" than SLy4. Varying 
the density, we found that the poles of the response func- 
tion at a;=0 occur only for po ^ 0.22 fm~^. This result 
does not guarantee that other types of instabilities could 
not develop |88| : however, we can rule out the most com- 
mon ones. 



3. Spherical shell structure 

The essence of the nuclear DFT is to be a global the- 
ory, whereby one unique functional (or family thereof) 
should be used to compute with reasonable accuracy var- 
ious properties of atomic nuclei from the lightest to the 
heaviest. Many of these properties depend on the single- 
particle shell structure. Figures S] and [5] show, respec- 
tively, the neutron s.p. energies in ^^Ca and proton s.p. 
energies in ^°^Pb obtained with UNEDFpre. They are 
compared with levels extracted from experiment [93| and 
those calculated with SLy4. In ^°®Pb. the overall agree- 
ment of the proton spectrum is very good. Furthermore, 
the neutron s.p. levels in ^°®Pb and proton and neutron 
levels in ^'^^Sn (not shown) agree well with experiment. 
As seen in Table |Vl the optimization produces a func- 
tional with an effective mass close to 1, which is proba- 
bly the reason the level density in ^'^^Pb and ^'^^Sn is well 
reproduced. Although the overall agreement for s.p. en- 
ergies is good, the systematic effect of high-j states being 
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FIG. 4: (color online) Neutron single-particle energies in Ca 
obtained from a HF calculation with the functional UNEDFpre. 
Experimental s.p. levels [931 ''■nd SLy4 results are shown for 
comparison. 
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slightly too high in energy is seen |57| . 

The neutron single-particle spectrum in ^®Ca is, how- 
ever, poorly reproduced. One of the most alarming fea- 
tures is the absence of the magic gap at iV=28 result- 
ing from a large s.p. level density and a reduced spin- 
orbit splitting. The s.p. proton spectrum of ^*Ca is only 
marginally better with the magic gap at Z=20 being too 
low, and the situation is similar in *'^Ca. 

The lack of observables directly probing s.p. properties 
(such as spin-orbit splittings or shell-gap sizes) in our 
objective function and the bias on heavy nuclei in the set 
of fit observables are undoubtedly the main reasons for 
the poor performance of UNEDFpre regarding the shell 
structure of light nuclei. Nevertheless, one must bear 
in mind that even when the optimization is exclusively 
focused on s.p. properties, standard Skyrme functionals 
perform poorly 



pro i 
&■ 



4-. Deformation properties 

The spherical shell structure determines many features 
of deformed nuclei. Indeed, the appearance of deformed 
states and shape coexistence effects can be related to s.p. 
levels and their couplings through symmetry-violating 
moments [93, [93| . Since the shell structure of light nuclei 
with UNEDFpre shows large deviations from experiment, 
it is interesting to test whether the new parameterization 
can nonetheless produce sensible deformation properties 
for medium-mass nuclei. To this end, we performed a se- 
ries of constrained HFB calculations for the sequence of 
Zr isotopes known to exhibit dramatic shape variations as 
a function of N . While nuclei near magic ^°Zr are known 
to be spherical, neutron- rich Zr isotopes with ^>100 pos- 
sess large prolate ground-state deformations, and 96-98 2j. 
exhibit a complex coexistence pattern [93, |9a] ■ On the 
proton-rich side, there is strong experimental evidence 
for large prolate deformation in iV=Z=40 system ^°Zr 

Figure [Hi shows the evolution of HFB-f LN deformation 
energy in the selected even-even Zr isotopes as a func- 
tion of the quadrupole deformation P2 ■ Each point was 
computed by imposing a constraint on the quadrupole 
moment {Q2) oc /32. Overall, the energy balance between 
spherical and deformed configurations is consistent with 
experiment. In particular, Zr isotopes with iV > 58 are 
correctly predicted to have prolate ground states coexist- 
ing with a secondary oblate minimum. Also, it is encour- 
aging to see that the ground state of ^'^Zr is predicted to 
be prolate, a feature that is not present in many Skyrme 
parameterizations [95|. 



Global mass table 



FIG. 5: (color online) Similar as in Fig. [J] but for proton 
single-particle energies in ^"^Pb. 



A good test of any EDF parameterization is its abil- 
ity to reproduce masses across the nuclear chart. Since 
our objective function contains the binding energies of a 
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FIG. 6: (color online) Deformation energy curves as functions 
of the quadrupole deformation P2 for selected even-even Zr 
isotopes calculated in the HFB+LN approach with UNEDFpre 
Skyrme functional. 



large set of nuclei, we expect good agreement with ex- 
perimental data, especially for heavy deformed systems. 

All even-even nuclei with N, Z > 8 have been calcu- 
lated with our two parameterizations according to the 
method presented in Q. Results have been posted 
for visualization and comparison with other EDF pa- 
rameterizations at http://massexplorer.org. The 
difference between experimental and theoretical bind- 
ing energies for the 520 even-even nuclei is shown in 
Fig. [71 An arclike trend [13 is seen for the SLy4 EDF; 
it has been attributed to an overemphasis on doubly 
magic nuclei during optimization. By contrast, both 
UNEDFpre and UNEDFnb show a much flatter behav- 
ior, while simultaneously reducing the mass residuals: 
RMSD(mass)==4.80MeV for SLy4, and 1.45 MeV and 
1.61 McV for UNEDFpre and UNEDFnb, respectively. 

To put things in perspective, we note that the best 
overall agreement with experimental masses obtained 
with the Skyrme EDF (on a larger data set that includes 
Ught and odd nuclei) is currently 0.582 MeV ^. How- 
ever, this excellent result was obtained at a price of sev- 
eral corrections on top of the EDF itself. In fact, a linear 
least-squares refit of the standard Skyrme EDF (also us- 
ing SLy4 as a starting point) to all even-even nuclear 
masses achieves a RMSD of around 1.7 MeV [Tol- Note 
also, that the RMSD for the masses of the UNEDFnb is 
higher by a 0.16 McV than for the UNEDFpre despite the 
larger domain available for parameter variation, which 
is due to the restricted set of masses used in this work. 
These figures suggest that UNEDFpre is probably within a 
few hundreds of keV of a globally optimal mass fit within 
the parameter space employed here. 

Close examination of Fig. [7] reveals that, while the 
global trend of binding energy errors has been improved, 
significant variations around that global trend still re- 
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FIG. 7: (color online) Binding energy residuals between the- 
ory and experiment for 520 even-even nuclei. The HFB-f LN 
results with SLy4 (top) are compared with those of UNEDFpre 
(middle) and UNEDFnb (bottom). 



main. To quantify this, we plot in Fig. |5] two-neutron 
separation energy residuals as a function of neutron num- 
ber N for the 520 nuclei of the previous set. The values 
of RMSD(S'2n) for SLy4 and UNEDFpre are, respectively, 
0.99 and 0.76 MeV, which indicate a significant improve- 
ment. If the set of 520 nuclei considered is divided into 
light {A < 80) and heavy (A > 80) subsets, the respec- 
tive RMSD(5'2„) values for SLy4 and UNEDFpre are 1.41 
and 1.45 MeV for light nuclei, and 0.85 and 0.45 MeV for 
heavy nuclei. This result stems from the bias toward 
heavy nuclei in our data set. 



6. Constraints on pairing strength from optimization 

Adjusting pairing interaction strengths represents a 
situation in which, by sequentially releasing a constraint 
on the EDF, one can dramatically improve the agree- 
ment with a subset of fit observablcs. The case in point 
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FIG. 8: (color online) Two-neutron separation energy residu- 
als between theory and experiment for 520 even-even nuclei. 
The HFB-I-LN results with SLy4 (top) are compared with 
those of UNEDFpre (middle) and UNEDFnb (bottom). 



is the interplay between pairing and shell structure. Since 
the shell correction to the binding energy favors low s.p. 
level density around the Fermi level, and the opposite 
is true for pairing contributions, an anticorrelation be- 
tween these two effects exists that results in a cancella- 
tion between shell and pairing energies [93| • If only total 
binding energies are subject to optimization, a reason- 
able fit can be obtained by, for example, an unphysical 
increase in pairing and a simultaneous unphysical varia- 
tion of s.p. shell structure. Indeed, since no data in our 
experimental data set directly probe s.p. energies, the 
lack of constraints on shell structure can dramatically 
impact pairing properties. 

Figure ini displays the OES residuals for three variants 
of calculations. In the first variant, the proton and neu- 
tron pairing strengths were kept equal and fixed at the 
standard value for SLy4 that yields an average neutron 
pairing gap in 
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FIG. 9: (color online) Neutron (left) and proton (right) OES 
residuals Ath — AeJp for the nuclei listed in Table IIVI The 
results with fixed (non-optimized) values of Vq" — Vq are 
marked by upside-down triangles. The optimized results are 
marked by triangles {V^ = Vjf ), dots (Vq" / V(f ; UNEDFnb), 
and squares (Vq" / Vq; UNEDFpre). 



^^°Sn equal to the experimental value of 



1.245 MeV [70|,|71|. In this case, the optimization proce- 
dure yields shell structure that resulted in overestimated 
pairing correlations, and the calculated RMSD for the 
OEM is 172 kcV. 

In the next step, we assumed proton and neutron pair- 
ing strengths to be identical V^ = Vq = Vq, and the 
constant Vq was included in the optimization set. The 
improvement on pairing energy was immediate, with the 
rms error on OES dropping down to 67keV. However, 
Fig. [5] clearly shows that OES for protons is almost sys- 
tematically underestimated. This observation calls for 
using different pairing strengths for neutrons and pro- 
tons, as was suggested by a recent large-scale survey |65l |. 

Our final optimization run was therefore carried out by 
considering independent strengths V^" and Vq in the fit. 
The rms error on OES has been further reduced to 57 keV 
in UNEDFnb and 59 kcV in UNEDFpre, and the two pairing 
strengths turn out to be significantly different; see Tables 
IVland rvTl Apart from possible global physics arguments, 
this result indicates that this optimization problem ben- 
efits from proton and neutron pairing strength being in- 
dependent parameters. 

We conclude this discussion with a word of warning: 
strictly speaking, the calculation of the OES requires 
computation of differences of binding energies. In odd 
nuclei, timc-rcvcrsal symmetry is broken, time-odd fields 
arc nonzero, and the ground-state should be computed as 
the lowest quasi-particlc excitation of a fully paired vac- 
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uum (blocking). Since the correct blocked state is not 
known beforehand, such calculations are much more in- 
volved than in even-even nuclei [7| . For this work, where 
the focus is on the optimization of the Skyrme functional 
itself and the pairing functional is limited, the extra cost 
of the proper treatment of odd nuclei was not deemed 
worth pursuing. 



B. Statistical Analysis of Optimization Results 

From a statistical viewpoint, our optimization prob- 
lem is also a nonlinear regression problem. For the true 
(but unknown) parameter value x* we define the errors 
between the theoretical value for the observable of type 
i in the nucleus j and its experimental counterpart as 



s,,j(x* 



(45) 



We assume every error e^ j is a random variable with 
expectation and that all Sij are independent and follow 
the same distribution. The optimization presented in 
Sec. IIIIAI estimates x, by the least-squares estimator 



argmin<j/(x) = -|lF(x)f 



(46) 



In the statistical setting, however, the random errors 
£ = {si_j} prevent the random variable x from always 
equaling x,. 



Confidence intervals 



Using the same notation as in (|4ip , we use a first-order 
approximation of the covariance matrix 



V^X^(x)K]g,g,^ 



Cov(x), (48) 



where parameters {gi}i=i....,nd ^^^ found by calculating 
central differences on the 2nx points {x ± 'r]kEk}k=i n ' 
where T]k > is chosen to be small. Although other ap- 
proximatio ns to the covariance matrix are possible, the 
authors of |lOl| state that V is their preferred approxi- 
mation because it is "simpler, less expensive, and more 
numerically stable" than alternative choices. 



TABLE VII: Optimal parameter values of UNEDFnb (no 
bounds), 95% confidence intervals, percentage of the initial 
guess for the scaling interval and standard deviation a. 



k 


Par. 


X 


95% CI 


% of Int. 


a 


1. 


Pc 


0.151046 


[0.149,0.153] 


10 


0.001 


2. 


^NM/^ 


-16.0632 


[-16.114,-16.013] 


5 


0.039 


3. 


^.NM 


337.878 


[302.692,373.064] 


70 


26.842 


4. 


NM 
^sym 


32.455 


[28.839,36.071] 


72 


2.759 


5. 


rNM 
-tJsym 


70.2185 


[11.108,129.329] 


296 


45.093 


6. 


1/m; 


0.95728 


[0.832,1.083] 


21 


0.096 


7. 




-49.5135 


[-55.786,-43.241] 


21 


4.785 


8. 


^pAp 


33.5289 


[-2.246,69.304] 


36 


27.292 


9. 


v^ 


-176.796 


[-194.686,-158.906] 


18 


13.648 


10. 


V,' 


-203.255 


[-217.477,-189.033] 


14 


10.850 


11. 


cr-' 


-78.4564 


[-85.137,-71.775] 


19 


5.097 


12. 


^PVJ 


63.9931 


[23.460,104.526] 


54 


30.921 



To characterize how the parameters change in a neigh- 
borhood of X* and x, we consider approximate confidence 
intervals. A 1 — a confidence interval flk C R is one in 
which we expect the true value Xk,* to lie 100(1 — a)% of 
the time, that is, with probability P(xk,* G ^k) = 1 — a. 

We note that the assumption of normally distributed 
residuals, e ~ A^(0, cr^/„^), made below, is the strongest 



one of this regression analysis. As pointed out in [2l| . 
theoretical (systematic) errors coming from an imperfect 
model are neither random nor generally independent, and 
their distribution is not rigorously normal. We carry out 
a standard analysis nonetheless in order to investigate 
constraints applied on our model. Therefore, the confi- 
dence intervals given here are to be understood as ranges 
of acceptable values for building parametcrizations of this 
particular model. 

Given normally distributed residuals and appropriate 
regularity conditions (as in [93], pages 23-25), a 1 — a 
confidence interval (CI) centered about Xk is 



Xk e 



\xk -ifel < yCov(x)fc,fc t 



rid—rixA- 



, (47) 



where i„^_„ 
distribution 



« is the 
with Ud - 



the covariance matrix Cov(x) 



- ^ quantilc of the t- 
dcgrces of freedom, and 

= E[(i-Ei)(i-Ei)^]. 



TABLE VIII: The same as Table |VTT1 except for the UN- 
EDFpre. 



k 


Par. 


X 


95% CI 


% of Int. 


O" 


1. 


Pc 


0.160526 


[0.160,0.161] 


10 


0.001 


2. 


^NM/^ 


-16.0559 


[-16.146,-15.965] 


45 


0.055 


3. 


^NM 


230 


- 


- 


- 


4. 


NM 
^sym 


30.5429 


[25.513,35.573] 


126 


3.058 


5. 


rNM 

-^sym 


45.0804 


[-20.766,110.927] 


219 


40.037 


6. 


1/m; 


0.9 


- 


~ 


- 


7. 


QpAp 


-55.2606 


[-58.051,-52.470] 


9 


1.697 


8. 


rjP^P 


-55.6226 


[-149.309,38.064] 


94 


56.965 


9. 


Vo" 


-170.374 


[-173.836,-166.913] 


3 


2.105 


10. 


vs 


-199.202 


[-204.713,-193.692] 


6 


3.351 


11. 


f^pVJ 


-79.5308 


[-85.160,-73.901] 


16 


3.423 


12. 


^PVJ 


45.6302 


[-2.821,94.081] 


65 


29.460 



Table I VIII shows the 95% confidence intervals, and 
standard deviations obtained when i~\k is chosen to be 
10~^ times the size of the scaling interval of parameter 
Xk- Standard deviations a are square roots of the diago- 
nal components of the covariance matrix Cov(x) and are 
often also referred to as errors of parameters. 
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Confidence intervals can therefore be valuable for test- 
ing the completeness of a given data set. In our case no 
data on giant resonances were included, which may ex- 
plain why a^J^ and L^J^^ remain imprecise. Similarly, our 
data set does not contain sufficiently many neutron-rich 
nuclei and/or entire isotopic sequences to pin down the 
isovector coupling constants. We also remark that the 
analysis based on confidence intervals is straightforward 
to perform once the (computationally-intensive) covari- 
ance matrix is known. 



2. Sensitivity Analysis 

The covariance matrix depends on the scaling of the 
parameters; hence, we will work with the standard cor- 
relation coefficient, 



^^^ Mass 



Proton radius 



OES 



Rk,l — 



Cov(a;fc,a;;) 
i/Var(xfe)Var(a;;) 



(49) 



which captures the (positive or negative) correlation be- 
tween parameters Xk and xi. Tables HXllXl provide the ap- 
proximate rij. X n^ correlation matrix R calculated when 
rjk is chosen to be 10~^ the size of the interval of in- 
terest of parameter Xk, for the solutions UNEDFpre and 
UNEDFnb, respectively. 

Overall, Tables lIXIIXl show that most parameters are 
interdependent, although the number of significant cor- 
relations with \Rkk\ > 0.8 is small. For UNEDFnb, where 
all parameters are free, we note that two pairs of NM 
parameters are well correlated: i^^M jg gy^ correlated 
withoc I56l.fl03 . while aT^, is 97% correlated with L?,m 
[53l 154 Il03j . The value of the (inverse of the) effective 
mass appears well correlated with both pairing strengths. 
We also notice strong correlation between the pairing 
strengths and the isoscalar spin-orbit coupling constants. 
Both observations reflect the interplay between single- 
particle level density and pairing discussed in Sec. IIV A 61 
We also notice that the proton pairing strength is signif- 
icantly correlated with the neutron pairing strength. 

In the case of the UNEDFpre parameterization, K^^ 
and 1/M* are removed from the sensitivity analysis. 
Nevertheless, we note that the various correlations be- 
tween parameters overall remain, even if they are atten- 
uated compared to the no-bound case. 

Next, we illustrate how sensitive the parameters Xk 
are to the different data types entering x^- masses, pro- 
ton radii, and OES. Here we focus only on the UN- 
EDFpre parameterization, as it is more realistic. We de- 
fine the rix X Hd Jacobian matrix J(x) of the residuals 



as J = (gi,...,g„ 



that is, the matrix formed as the 



juxtaposition of the n^ column- vectors gi {nd being the 
number of data, n^ the number of parameters). The 
nx X nd sensitivity matrix S is 



5(x)=[j(x)J^(x)]"'j(x) 



(50) 



For each line in the sensitivity matrix (each parame- 
ter), we can compute the partial sums over each of the 
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FIG. 10: (color online) Sensitivity of the parameters of UN- 
EDFpre to different data types entering x^- The EDF param- 
eters are labeled as in Table IVIII 



three types of data. This coinputation gives us a measure 
of the change of the parameter under a global change of 
all the data of a given type. Figure [TOl shows the relative 
change of parameter xu when such an average datum of 
an observable is changed. For example, for i= 1 (masses), 
it shows the change in Xk under a variation of all exper- 
imental masses. 

All of the bars in Fig. [10] have been renormalized to 
unity, and only relative strengths between mass, radii, 
and OES data are shown. A large percentage contribu- 
tion from data type i means that Xk is very sensitive to 
changes in i, and other data types have little impact on it 
at the convergence point. As expected, pairing strengths 
(parameters 9 and 10) are primarily affected by OES 
data. It is worth noting the very similar sensitivity of 
the spin-orbit coupling constants (parameters 11 and 12) 
on all 3 types of data. Also, nuclear matter parameters 
appear to be significantly more dependent on the proton 
radius than other coupling constants. This is not sur- 
prising, considering the relation between the saturation 
density and the WignerSeitz radius. 

The integrated information contained in Fig. llOl cannot 
assess the impact of a particular data piece on model pa- 
rameters; hence, a more detailed analysis is needed. To 
this end, for each experimental observable dij (masses, 
radii, OES), we compute the global change in x as in- 
dividual data di^j change (one-by-one) by Clm^, namely, 
200 keV for masses, 0.002 fm for proton radii, and 50 
keV for OES. In this way we can, for instance, evalu- 
ate the possible importance of some new experimental 
observable on a given model [3a| . Figure [TT] shows the 
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TABLE IX: Correlation matrix (|49|) for the UNEDFnb parameter set (no bounds). 



Pc 


1.00 
























E^^/A 


-0.04 


1.00 






















^NM 


-0.87 


0.16 


1.00 




















NM 
"'sym 


-0.05 


-0.72 


-0.29 


1.00 


















rNM 
-^sym 


-0.09 


-0.62 


-0.23 


0.97 


1.00 
















l/Mt 


-0.05 


0.05 


0.09 


-0.10 


-0.10 


1.00 














nP^p 


-0.23 


0.24 


0.34 


-0.25 


-0.20 


-0.86 


1.00 












QpAp 


-0.22 


0.29 


0.34 


-0.65 


-0.76 


-0.08 


0.28 


1.00 










Vo" 


0.02 


-0.02 


-0.06 


0.06 


0.06 


-0.99 


0.87 


0.12 


1.00 








V,' 


0.01 


-0.14 


-0.10 


0.26 


0.27 


-0.95 


0.78 


-0.07 


0.93 


1.00 






C/?^^ 


0.07 


-0.03 


0.04 


-0.14 


-0.17 


-0.72 


0.78 


0.32 


0.73 


0.65 


1.00 




^PVJ 


-0.07 


-0.35 


-0.12 


0.58 


0.66 


0.06 


-0.26 


-0.64 


-0.08 


0.05 


-0.38 


1.00 




Pc 


^NM/^ 


j^^M 


NM 
i^syni 


rNM 


1/AC 


nP^p 

^0 


nP'^p 




VS 


^pWJ 


^pyj 



TABLE X: Correlation matrix (|49|) for the UNEDFpre parameter set. 



Pc 


1.00 
























^NM/^ 


-0.28 


1.00 






















7V-NM 


- 


- 


- 




















NM 
'^sym 


-0.10 


-0.88 


- 


1.00 


















rNM 
-'^sym 


-0.17 


-0.80 


- 


0.97 


1.00 
















1/M* 


- 


- 


- 


- 


- 


- 














r^pi^p 

^0 


0.09 


0.80 


- 


-0.81 


-0.74 


- 


1.00 












fjpAp 


0.20 


0.35 


- 


-0.47 


-0.66 


- 


0.23 


1.00 










Vo" 


0.02 


0.21 


- 


-0.23 


-0.25 


- 


0.23 


0.23 


1.00 








v,f 


-0.13 


-0.42 


- 


0.52 


0.56 


- 


-0.29 


-0.45 


-0.14 


1.00 






cC' 


0.37 


-0.14 


- 


0.02 


-0.00 


- 


0.44 


-0.02 


0.09 


0.16 


1.00 




cf' 


-0.06 


-0.18 


- 


0.27 


0.33 


- 


-0.38 


-0.20 


-0.01 


0.00 


-0.37 


1.00 




Pc 


^NM/^ 


^.NM 


NM 

^sym 


rNM 

-'^sym 


l/Mt 


nP'^p 


j^pAp 


VS 


VE 


ct^-' 


ct' 



quantity 



(5a;/cr|| = 



\ 






(51) 



with (5xfc being the change in the value of the parameter 
Xfc under a change of the data d^.j, for all 71^=108 data 
points. This is nothing but the norm of the total change 
in units of the standard deviation a}^ , defined as before by 
(Tfc = \/Cov(a;fe,xs;). Large changes in x mean that the 
parameter values are highly sensitive to the particular 
value of di^y 

In principle, the sensitivity analysis can be performed 
at any point x of the ri3;-dimensional parameter space: 
for a given scalar function /(x) of the type (pS]) . the 
sensitivity at point x is only based on the local gradient. 
In particular, it is totally independent of the procedure 
that leads to the specific selection of x. It only depends 
on the degrees of freedom, i.e., free parameters, retained 
in/(x). 

The question of degrees of freedom is highly relevant 
in the context of the UNEDFpre parameter set, where 
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two parameters are actively constrained at the solution. 
These constraints have been directly implemented by re- 
stricting the domain where the y^ function is evaluated, 
and not by modifying the function by adding a penalty. 

The net result of imposing the bound constraints is 
that only 10 parameters out of twelve are allowed to 
change near the end of the optimization process. One 
has, therefore, two options as far as the sensitivity anal- 
ysis is concerned: (i) Remove these two parameters from 
the set of active parameters and calculate the Jacobian 
J(x) and the sensitivity matrix S'(x) with only rir^ — 2=10 
parameters; (ii) Keep all n^=Vl parameters in the calcu- 
lation of the Jacobian and doing a tangent plan approx- 
imation to obtain the relevant covariance and sensitivity 
matrices. In this way, the uncertainties of the other 10 
parameters are affected by the local fluctuations of the 
surface induced by these two actively constrained param- 
eters. 

The alternative (i) boils down to computing the gradi- 
ent at point X in a 10-parameter subspace of the original 
12-parameter space. In this subspace, x is the stable 
point of the x^ function. The curve labeled UNEDFpreio 
in Fig. [TT] corresponds to this approach. This sensitiv- 
ity response is compared with that performed at the free 
minimum x''^'^^ of ^ in the full 12-parameter space of 
UNEDFnb. Since both sets correspond to (unconstrained) 
minima in their respective spaces, the overall changes 
in X are of the same order of magnitude and they are 
very small, ||(5a;/(T|| ~ 0.01. This indicates that the set 
of fit observables has been chosen very consistently. In- 
deed the mass of deformed ^'^^Fm is the single observable 
that yields the noticeable parameter variations around 
the UNEDFnb minimum while in the case of UNEDFprcio 
no sensitivity to a single piece of data can be noticed. 

By contrast, within the alternative (ii), x is not an 
unconstrained minimizer in the full 12-parameter space. 
In such an approach, we find the overall sensitivity to be 
about 2 orders of magnitude larger than what is depicted 
in Fig. [TT] This reflects the fact that 2:3 and xg are far 
away from the unconstrained minimum in that space, at 
the same time being strongly correlated with other pa- 
rameters. In this case, the sensitivities depend strongly 
on the actual value of x and the way the domain is con- 
strained. For this reason, the option (ii) is of no practical 
interest in the comparison with the no-bounds results. 



V. CONCLUSIONS 

One of the major challenges for the low-energy nuclear 
theory is to construct the global nuclear energy density 
functional of spectroscopic quality, rooted in microscopic 
theory. An important element of this program is to opti- 
mize the parameters of the functional on a set of experi- 
mental observables and selected theoretical pseudo-data. 
This work shows how such an optimization can be done 
by using modern optimization algorithms and nonlinear 
regression analysis. 



The purpose of this study was to optimize the stan- 
dard Skyrme functional based on a set of experimen- 
tal data (masses, charge radii, and odd-even mass dif- 
ferences) pertaining to 72 spherical and deformed nuclei 
amenable to a mean-field description. The new model- 
based optimization algorithm POUNDerS was compared 
with other standard derivative-free optimization meth- 
ods such as Neldcr-Mead and was found to be signifi- 
cantly better in terms of speed, accuracy, and precision. 

The optimization was carried out at the fully 
self-consistent, deformed Hartree-Fock-Bogoliubov level. 
Here, we took advantage of the efficient DFT solver HF- 
BTHO optimized in the first phase of the project. We 
have implemented various improvements that enable us 
to quickly compute global self-consistent mass tables. 
This capability is essential for optimization. 

As a result of the twelve-parameter optimization of 
Skyrme EDF, we arrived at two solutions. The first 
one corresponds to a minimum (stable point) in the con- 
sidered parameter space. The corresponding functional 
UNEDFnb describes well the assumed set of fit observ- 
ables, but its incompressibility parameter is too large, as 
this property has not been well constrained by our data 
set. The second optimization was carried out assum- 
ing hard bounds on the nuclear matter parameters. For 
the bound-constrained solution, the nuclear incompress- 
ibility and scalar effective mass appear at their respec- 
tive bounds. The resulting parameter set UNEDFpre gives 
good agreement with experimental masses, radii, and de- 
formations and seems to be free of finite-size instabilities. 
In particular, for two-neutron separation energies and 
masses of even-even heavy nuclei with A > 80, UNEDFpre 
yields the rms deviation of 0.45 MeV and 1.2 MeV, re- 
spectively, which is a satisfying result. We emphasize 
that the original Skyrme EDFs seem to be inherently 
limited in this respect, as demonstrated in [lO[, unless 
specific corrections are introduced. Our result is there- 
fore in line with the best expectations one could have 
for such EDFs. Nevertheless, the lack of specific con- 
straints on the shell structure in our data set implies 
that single-particle levels of light nuclei are not well re- 
produced. For that reason, UNEDFpre may not yet be rec- 
ommended for truly global applications across the chart 
of the nuclides. However, this functional is expected to 
work well for heavy nuclei and should be considered as 
a reference against which more advanced EDFs will be 
benchmarked. 

We have also applied full-fiedged regression diagnostics 
on UNEDFnb and UNEDFpre, focusing on statistical cor- 
relations between ED parameters and the sensitivity of 
parameters to variations in fit observables. To this end, 
we computed and analyzed the correlation and sensitivity 
matrices at the optimal parameter set. This kind of non- 
linear regression analysis is expected to be helpful when 
designing next generation EDFs. Moreover, the statisti- 
cal tools presented in this study can be used to pinpoint 
specific nuclear observables that are expected to strongly 
affect the developments of the nuclear universal density 
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functional. 
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